//! A solid convex polygon. It is assumed that the interior of the polygon is to
//! the left of each edge.

const std = @import("std");
const Vec2 = @import("../vec2.zig");
const Point = @import("../point.zig");
const Common = @import("common.zig");
const AABB = Common.AABB;
const MassData = Common.MassData;
const Proxy = @import("proxy.zig");
const Hull = @import("hull.zig");

pub const Polygon = @This();

/// The polygon vertices
vertices: []Point,

/// The outward normal vectors of the polygon sides
normals: []Vec2,

/// The centroid of the polygon
centroid: Point,

/// The external radius for rounded polygons
radius: f32,

allocator: std.mem.Allocator,

/// Offset Polygon: radius=0 and transform=null
/// Offset Rounded Polygon: radius>0 and transform != null
pub fn init(hull: *const Hull, radius: f32, tx: ?*const Point.Transform, allocator: std.mem.Allocator) *Polygon {
    var self = allocator.create(Polygon) catch @panic("Out of memory");
    self.allocator = allocator;
    self.radius = radius;

    self.vertices = allocator.alloc(Point, hull.points.len) catch @panic("Out of memory");
    if (null == tx) {
        std.mem.copyForwards(Point, self.vertices, hull.points);
    } else {
        for (0..hull.points.len) |i| {
            self.vertices[i] = hull.points[i].ifTransform(tx.?);
        }
    }

    self.normals = allocator.alloc(Vec2, hull.points.len) catch @panic("Out of memory");
    for (0..hull.points.len) |i| {
        const i_next = if (i + 1 < hull.points.len) i + 1 else 0;
        const edge = self.vertices[i].vectorTo(&self.vertices[i_next]);
        self.normals[i] = edge.crossVS(1.0).normalize();
    }

    self.centroid = computeCentroid(self.vertices);

    return self;
}

pub fn deinit(self: *Polygon) void {
    self.allocator.free(self.vertices);
    self.allocator.free(self.normals);

    self.allocator.destroy(self);
}

pub fn makeProxy(self: *const Polygon) *Proxy {
    return Proxy.init(self.vertices, self.radius, self.allocator);
}

pub fn isInnerPoint(p: *const Point) bool {
    _ = p;
}

pub fn transform(self: *Polygon, tx: *const Point.Transform) void {
    for (0..self.vertices.len) |i| {
        self.vertices[i].transform(tx);
        self.normals[i].rotate(tx.rotation);
    }
    self.centroid.transform(tx);
}

pub fn ifTransform(self: *const Polygon, tx: *const Point.Transform) *Polygon {
    var new_one = self.allocator.create(Polygon) catch @panic("Out of memory");

    new_one.allocator = self.allocator;
    new_one.radius = self.radius;

    new_one.vertices = self.allocator.alloc(Point, self.vertices.len) catch @panic("Out of memory");
    new_one.normals = self.allocator.alloc(Vec2, self.vertices.len) catch @panic("Out of memory");
    for (0..self.vertices.len) |i| {
        new_one.vertices[i] = self.vertices[i].ifTransform(tx);
        new_one.normals[i] = self.normals[i].ifRotate(tx.rotation);
    }
    new_one.centroid = self.centroid.ifTransform(tx);

    return new_one;
}

/// Polygon mass, centroid, and inertia.
/// Let rho be the polygon density in mass per unit area.
/// Then:
/// mass = rho * int(dA)
/// centroid.x = (1/mass) * rho * int(x * dA)
/// centroid.y = (1/mass) * rho * int(y * dA)
/// I = rho * int((x*x + y*y) * dA)
///
/// We can compute these integrals by summing all the integrals
/// for each triangle of the polygon. To evaluate the integral
/// for a single triangle, we make a change of variables to
/// the (u,v) coordinates of the triangle:
/// x = x0 + e1x * u + e2x * v
/// y = y0 + e1y * u + e2y * v
/// where 0 <= u && 0 <= v && u + v <= 1.
///
/// We integrate u from [0,1-v] and then v from [0,1].
/// We also need to use the Jacobian of the transformation:
/// D = cross(e1, e2)
///
/// Simplification: triangle centroid = (1/3) * (p1 + p2 + p3)
///
/// The rest of the derivation is handled by computer algebra.
pub fn computeMass(self: *const Polygon, density: f32) MassData {
    const count = self.vertices.len;
    const radius = self.radius;

    var vertices = self.allocator.alloc(Point, count) catch @panic("Out of memory");
    defer self.allocator.free(vertices);

    if (radius > 0.0) {
        // Approximate mass of rounded polygons by pushing out the vertices.
        const sqrt2_radius: f32 = 1.412 * radius;

        for (0..count) |i| {
            const j = if (0 == i) count - 1 else i - 1;

            var mid = self.normals[i].ifAdd(&self.normals[j]).normalize();
            mid.mul(sqrt2_radius);

            vertices[i] = self.vertices[i].ifMove(&mid);
        }
    } else {
        std.mem.copyForwards(Point, vertices, self.vertices);
    }

    var center = Point.init(0, 0);
    var area: f32 = 0.0;
    var rotational_inertia: f32 = 0.0;

    // Get a reference point for forming triangles.
    // Use the first vertex to reduce round-off errors.
    const r = vertices[0];

    const inv3: f32 = 1.0 / 3.0;

    for (1..count - 1) |i| {
        const p = vertices[i];
        const p_next = vertices[i + 1];

        const e1 = r.vectorTo(&p);
        const e2 = r.vectorTo(&p_next);

        const D: f32 = e1.cross(&e2);

        const triangle_area: f32 = 0.5 * D;
        area += triangle_area;

        // Area weighted centroid, r at origin
        center.move(&e1.ifAdd(&e2).ifMul(triangle_area * inv3));

        const e1x = e1.x;
        const e1y = e1.y;
        const e2x = e2.x;
        const e2y = e2.y;

        const intx2: f32 = e1x * e1x + e2x * e1x + e2x * e2x;
        const inty2: f32 = e1y * e1y + e2y * e1y + e2y * e2y;

        rotational_inertia += (0.25 * inv3 * D) * (intx2 + inty2);
    }

    const mass: f32 = density * area;
    const inv_area: f32 = 1.0 / area;
    center.x *= inv_area;
    center.y *= inv_area;
    const mass_center = center.ifMoveAwayOrigin(r.x, r.y);

    // Inertia tensor relative to the local origin (point s).
    var _rotational_inertia = density * rotational_inertia;
    // Shift to center of mass then to original body origin.
    _rotational_inertia += mass * (mass_center.x * mass_center.x + mass_center.y * mass_center.y - (center.x * center.x + center.y * center.y));

    return MassData.init(mass, mass_center, _rotational_inertia);
}

pub fn computeAABB(self: *Polygon, tx: *const Point.Transform) AABB {
    var lower = self.vertices[0].ifTransform(tx);
    var upper = lower;

    for (1..self.vertices.len) |i| {
        const v = self.vertices[i].ifTransform(tx);
        lower = lower.newMinPoint(&v);
        upper = upper.newMaxPoint(&v);
    }

    lower = lower.ifMove2Origin(self.radius, self.radius);
    upper = upper.ifMoveAwayOrigin(self.radius, self.radius);

    return AABB.init(lower, upper);
}

/// 质心法
/// 1、将凸多边形划分为n−2个三角形(从某一顶点出发)
/// 2、计算每个三角形的面积Ai和质心(Cxi,Cyi)
/// 3、多边形质心为各三角形质心的加权平均
/// x = sum(Ai*Cxi)/sum(Ai), 其中 i in [1,n-2]
/// y = sum(Ai*Cyi)/sum(Ai), 其中 i in [1,n-2]
///
/// 公式可简化为：
/// A = 0.5 * sum(xi*y(i+1) - x(i+1)*yi), 其中 i in [0,n-1]
/// x = sum((xi+x(i+1)) - (xi*y(i+1) - x(i+1)*yi))/6A, 其中 i in [0,n-1]
/// y = sum((yi+y(i+1)) - (xi*y(i+1) - x(i+1)*yi))/6A, 其中 i in [0,n-1]
fn computeCentroid(vertices: []Point) Point {
    var x: f32 = 0.0;
    var y: f32 = 0.0;
    var area: f32 = 0.0;

    for (0..vertices.len) |i| {
        const p = vertices[i];
        const p_next = vertices[if (i + 1 < vertices.len) i + 1 else 0];

        const _area = p.x * p_next.y - p_next.x * p.y;
        area += _area;
        x += (p.x + p_next.x) * _area;
        y += (p.y + p_next.y) * _area;
    }
    area *= 3;

    return Point.init(x / area, y / area);
}
